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We study the structure, stability, and dynamics of dark solitary waves in parabolically trapped, 
collisionally inhomogeneous Bose-Einstein condensates with spatially periodic variations of the scat- 
tering length. This collisional inhomogeneity yields a nonlinear lattice, which we tune from a small- 
amplitude, approximately sinusoidal structure to a periodic sequence of densely-spaced spikes. We 
start by investigating time- independent inhomogeneities, and we subsequently examine the dy- 
namical response when one starts with a collisionally homogeneous BEC and then switches on an 
inhomogeneity either adiabatically or nonadiabatically. Using Bogoliubov-de Gennes linearization 
as well as direct numerical simulations of the Gross-Pitaevskii equation, we observe dark solitary 
waves, which can become unstable through oscillatory or exponential instabilities. We find a critical 
wavelength of the nonlinear lattice that is comparable to the healing length. Near this value, the 
fundamental eigenmode responsible for the stability of the dark solitary wave changes its direction 
of movement as a function of the strength of the nonlinearity. When it increases, it collides with 
other eigenmodes, leading to oscillatory instabilities; when it decreases, it collides with the origin 
and becomes imaginary, illustrating that the instability mechanism is fundamentally different in 
wide- well versus narrow- well lattices. When starting from a collisionally homogeneous setup and 
switching on inhomogeneities, we find that dark solitary waves are preserved generically for aligned 
lattices. We briefly examine the time scales for the onset of solitary- wave oscillations in this scenario. 



I. INTRODUCTION 

Because of advances in the field of Bose-Einstein con- 
densation, solitary waves have received considerable at- 
tention in atomic physics during the past 15 years jl|, Q. 
In particular, many studies have examined a mean- field 
description of atomic Bose-Einstein condensates (BECs) 
using the Gross-Pitaevskii (GP) equation. The GP equa- 
tion includes a cubic nonlinearity, with a strength pro- 
portional to the s-wave scattering length, which arises 
from the interatomic interactions. This mean- field de- 
scription has been important for theoretical and exper- 
iment investigations of phenomena such as bright soli- 
tary waves [Jlsl , dark solitary waves [6l-[l3j , gap solitary 
waves [HI, multi-component solitary waves fl5l-[2lj. and 
Faraday waves [22-24]. Solitary waves have been studied 
in great detail in the presence of various external poten- 
tials, including linear periodic lattices (so-called "optical 
lattices") and harmonic traps. This has yielded tremen- 
dous insights into a large variety of phenomena, includ- 
ing Bloch oscillations, Landau-Zener tunneling, modu- 
lational ( "dynamical" ) instabilities and gap excitations, 
and more jgfl, EH, \2^2^ . 

Numerous techniques have been developed that enable 
the control and manipulation of coherent states in BECs. 
These include static (homogeneous and inhomogeneous) 
electric and magnetic fields [29| , optical devices [30] , and 
near-field radio- frequency devices [31]. One can vary a 
BECs external (trapping) potential while independently 
and simultaneously changing the strength of the nonlin- 
earity by tuning interatomic interactions. The interac- 
tion among the atoms can be adjusted experimentally 



over a very broad range by employing either magnetic 
[32l |33( | or optical Feshbach resonances [34] , and the ma- 
nipulation of BECs using Feshbach resonances has led to 
numerous insights. Experimental achievements include 
the formation of bright solitary waves and solitary- wave 
trains for 7 Li 0, 3 and 85 Rb [5] atoms by tuning the 
interatomic interaction within a stable BEC from repul- 
sive to attractive, the formation of molecular condensates 
[35^ and the probing of the BEC-BCS crossover [36| . 
Theoretical studies include the prediction that a time- 
dependent modulation of the scattering length can be 
used to stabilize attractive two-dimensional (2D) BECs 
against collapse j37j or to create robust matter-wave 
breathers in ID BECs [38]. Temporal modulation of the 
GP equation's nonlinearity and its impact on collapse 
properties and on other nonlinear phenomena (such as 
modulational instabilities) have been examined system- 
atically in nonlinear optics [391 ]. 

Atomic matter waves also exhibit novel features un- 
der the influence of a spatially varying scattering length, 
which yields a spatially-dependent nonlinearity coeffi- 
cient in the GP equation. Numerous works have consid- 
ered matter- wave dynamics in such "collisionally inhomo- 
geneous" environments. Theoretical predictions include 
adiabatic compression of matter waves [H, Q •> enhance- 
ment of the transmitivity of matter waves through bar- 
riers |45|, |46| , dynamical trapping of solitary waves [45| , 
and a derealization transition of matter waves [I?} • Lin- 
ear [43|, HH, parabolic |48|, r andom [49| , periodic (i.e., 
nonlinear lattices) |47|, I5CH52] , and localized (step-like) 
[53l-[55j inhomogeneities have all been considered. There 
have also been several detailed mathematical studies [56l — 
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[58[ as well as examinations of analogous situations in op- 
tics [59j| . Additionally, the interplay between linear and 
nonlinear lattices has been examined in both continuum 
[60| and discrete [6l[ settings. Two important recent de- 
velopments on BECs in collisionally inhomogeneous en- 
vironments include a broad review of relevant theoretical 
activity [62| and the experimental creation of a spatially 
periodic effective nonlinearity in a Yb BEC [63|. 

In this paper, which is motivated in part by the re- 
cent experimental work of Ref. [63|, we consider spa- 
tially periodic scattering lengths that can be tuned us- 
ing Feshbach resonances from small-amplitude, approx- 
imately sinusoidal structures to a periodic sequence of 
densely-spaced spikes. We consider "aligned" and "anti- 
aligned" periodic structures and compare and contrast 
the solitary- wave dynamics for each case. We first con- 
sider situations in which collisional inhomogeneities are 
always present (i.e., in which the coefficient of the nonlin- 
earity is time-independent) before moving on to ones in 
which they are turned on either adiabatically or abruptly. 

When considering time-independent inhomogeneities, 
we employ both direct numerical simulations of the GP 
equation and a Bogoliubov-de Gennes (BdG) analysis 
to examine the linear stability of dark solitary waves. 
We thereby identify two distinct mechanisms of instabil- 
ity for the dark solitary waves — an exponential one and 
an oscillatory one — and we examine when each of these 
arises. Importantly, we find and quantify a change of 
instability mechanism as one switches between wide-well 
and narrow- well lattices. When switching on inhomo- 
geneities in time, we find for aligned lattices that the 
dark solitary waves are generically robust, but that addi- 
tional excitations in the form of gray solitary waves also 
emerge (especially when the inhomogeneities are turned 
on nonadiabatically ) . We briefly examine the time scales 
for the onset of solitary- wave oscillations in this scenario. 

The rest of our presentation is organized as follows. In 
Section [IIJ we present the GP equation and our setup of 
spatially periodic scattering lengths. We consider time- 
independent scattering lengths and discuss the dark soli- 
tary waves arising in this situation in Section lllll We then 
consider the dynamical response of switching on the col- 
lisional inhomogeneities in Section IIVI Finally, we con- 
clude in Section [V] and present some directions for future 
work. 



II. SETUP AND MODEL 

The dynamics of a cigar-shaped BEC can be approx- 
imated in the mean field by the quasi- ID GP equation 

mam 



and g(z,t) is a (suitably normalized) spatially and tem- 
porally modulated coefficient. We suppose that the BEC 
is in a harmonic trap, so V(z,i) = j^jZ 2 . In Eq. (pQ), we 
have rescaled the condensate density in units of Bhuj z /\q\, 
length in units of a z /\/2B, time in units of (Bujz) -1 , and 
energy in units of HBuj z . The parameter a z — \Jhj (muj z ) 
is the axial oscillator length, h is Planck's constant, m is 
the mass of an atom in the BEC, and uj z is the axial trap 
frequency. Additionally, q = 2aftw± is the dimensional 
nonlinearity coefficient, where a is the s- wave scattering 
length and we have averaged over the transverse direc- 
tions (assuming a ground-state wave profile) [20|, Hl| . We 
choose B = 10 for convenience. We consider a 87 Rb BEC 
as an example and consequently use the parameter values 
a w 5.5 x 10 _9 m, uj z — 2tt x 4s -1 , and uj± = S0uj z . 

In our calculations, the nondimensional numbers of 
atoms N = J^°° \u\ 2 dz is roughly 500. The number 
of atoms J\f = (h 3/2 [Bcj z /2m] 1/2 /\q\)N is thus roughly 
7000, though the phenomena that we observe are ro- 
bust over different values of N. The two criteria for 
the quasi- ID regime are satisfied [65]: (1) the trap is 
highly anisotropic, as Q = oj z /oj± — 0.0125 <C 1; and 
(2) the parameter d = JVQa/a± « 0.8 < 1, where 
a± = ^/H/(muj±) is the oscillator length in the radial 
direction. 

The normalized coefficient g oc a can be either posi- 
tive or negative. The sign of g depends on the atomic 
species; repulsive interatomic interactions yield g > 0, 
and attractive ones yield g < 0. The sign and mag- 
nitude of g can both be changed using Feshbach reso- 
nances, which make it possible (in principle) to manipu- 
late the sign and strength of atomic interactions [H, [33| . 
It is also possible to vary the scattering length and con- 
sequently the nonlinearity coefficient g in space and time 
by tuning an external field in the vicinity of a Feshbach 
resonance. As we noted previously, spatial variation of 
the 5-wave scattering length using Feshbach resonances 
was recently demonstrated experimentally [63]. Tempo- 
ral variations of BEC scattering lengths using Feshbach 
resonances have been studied for more than 10 years. 

In this paper, we consider repulsive BECs with a spa- 
tially periodic nonlinearity coefficient g = g(z,t). We 
assume that the nonlinearity does not change sign. We 
start by considering time-independent s-wave scattering 
lengths, for which g(z , t) = g(z). We consider the exper- 
imentally realistic [63| functional form 



Ag(z) = 



g + Ag(z) , 

gm 

1 + g s sin(/cz + 4>) 



(2) 



■ 9 . , 



a_ 2 

dz 



2 + V(z,t) +g(z,t)\u(z,t)\ 2 - n ) u(z,t) 



(1) 



where u(z,t) is the macroscopic wave function, /i is the 
chemical potential [67j, V(z,t) is an external potential, 



where the wavenumber k determines the wavelength of 
the nonlinear lattice g(z), the parameter go = 1 (by nor- 
>malization), and Ag(z) > 0. We let g m G [0,1], and we 
use g s G [0, 0.95] so that Ag(z) is always finite. In our nu- 
merical simulations, we consider wave numbers k G [0, 15]. 

The parameter g m determines the magnitude of the 
spatial modulation of g(z), and the parameter g s deter- 
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mines the shape of the lattice, which is approximately 
sinusoidal when g s <C 1 and resembles a periodic train 
of thin spikes as g s —> 1. We avoid a singularity in g(z) 
by bounding g s away from 1. The angle is a constant 
and determines a shift in alignment between the nonlin- 
ear lattice and the harmonic potential V(z). We consider 
two situations: (1) aligned nonlinear lattices, for which 
cj) = ^ 7 so the minimum of the harmonic potential V(z) 
(which is located at z = 0) coincides with a minimum of 
g(z)\ and (2) anti-aligned lattices, for which <ft = — ^, so 
the minimum of the harmonic trap coincides with a maxi- 
mum of g(z). As we discuss later, the stability properties 
of solitary waves differ in the two cases. 

In the sections below, we examine the existence, stabil- 
ity, and dynamics of dark solitary waves as the strength 
and shape of the nonlinear lattice g(z) is varied. We es- 
tablish existence by finding standing- wave solutions u(z), 
for which the right-hand side of Eq. (pQ) vanishes. We then 
examine the stability of these solutions via linearization 
around u. In other words, we perform a BdG analysis: 
we consider an 0(e) correction to the GP equation (pQ) 
by writing 

u(z, t) = u(z) + e [a(z)e- iujt + b*(z)e iujt ] , 

where the asterisk * denotes complex conjugation. We 
examine the dynamics of the standing wave by perturbing 
it and computing its temporal evolution using Eq. (pQ). 
We then study the response of solitary waves to both 
gradual (adiabatic) and abrupt (nonadiabatic) changes 
of the atomic interactions. Both types of changes can, in 
principle, be introduced experimentally using a Feshbach 
resonance. 



III. DARK SOLITARY WAVES 

As we will discuss in this section, the stability of the 
observed dark solitary waves depends on the parameter 
values of the nonlinear lattice. In Figs. [I] and [2j we show 
the dynamics for several situations in aligned and anti- 
aligned lattices. We show solitary- wave solutions of the 
stationary GP equation (obtained via Newton iteration) 
in the top panels, where we also display the nonlinear 
lattice and the harmonic trap. The dip in the wavefunc- 
tions near z = helps illustrate that these are in fact 
dark solitary waves. We show their corresponding eigen- 
frequencies (computed using the BdG equations) in the 
middle panels, and we show the temporal evolution of the 
solitary waves using direct numerical simulations of the 
time-dependent GP equation (pQ) in the bottom panels. 

We observed two qualitatively distinct situations for 
both aligned and anti-aligned nonlinear lattices: unsta- 
ble dark solitary waves with one pair of purely imaginary 
eigenfrequencies and unstable dark solitary waves with 
a quartet of complex eigenfrequencies. As the value of 
the wavenumber k increases, the solution in the aligned 
lattice undergoes a transition from windows of mild os- 
cillatory instabilities, which are indicated by quartets of 



complex eigenfrequencies that result from Hamiltonian- 
Hopf bifurcations, to strong instabilities, which are indi- 
cated by purely imaginary frequencies that result from 
crossing the origin of the spectral plane. Interestingly, 
we observe the opposite dependence on k in the anti- 
aligned lattice. We illustrate these small- k (left panels) 
and large- k (right panels) results for aligned lattices in 
Fig. [D and for anti-aligned lattices in Fig. [2j 

The dichotomy between the dynamics for aligned and 
anti-aligned lattices indicates that, given a fixed set of 
parameter values, one can in principle shift the nonlinear 
lattice (or, equivalently, the magnetic trap) and control 
the strength of the instability of a dark solitary wave. 
This makes it possible to transition between a regime of 
strong instability and a regime of alternating windows of 
weak instability. From the spatio-temporal evolution of 
the dark solitary waves, which we illustrate in the bot- 
tom panels, we observe in both cases that the waves ul- 
timately become displaced from the center and oscillate 
between two turning points in the parabolically trapped 
BEC. In the small- k regime, the strong interaction of the 
propagating dark solitary-wave structures leads to their 
eventual decay after only a few oscillations. When k is 
large, however, the BEC appears to reach a "homoge- 
nized" limit that is characterized by a reduced range of 
spatial variation (due to the nonlinear lattice) within the 
trap. Ultimately, as k increases further, the role of the 
nonlinear lattice is to yield an effective averaged modu- 
lation of the 5-wave scattering length. 

Given a regime of the wavenumber /c, we vary g m from 
by a small amount and monitor the change of eigen- 
frequency spectrum that we obtain using the BdG anal- 
ysis. In particular, we examine the four (pairs of) eigen- 
frequencies of smallest magnitude. We illustrate typi- 
cal trends of the eigenfrequency spectrum variation in 
Figs. [3H6l In order, these figures show results for the 
small- k regime for a BEC in an aligned nonlinear lat- 
tice, the large- k regime for an aligned lattice, the small- k 
regime for an anti-aligned lattice, and the large- k regime 
for an anti-aligned lattice. In all eigenfrequency com- 
putations, we are concerned predominantly with the so- 
called anomalous mode (i.e., the only mode that has neg- 
ative Krein signature [65]). The anomalous mode is lo- 
cated at u = ±l/(By/2) for g m = 0, and we plot it using 
red circles in the figures. In the small- k regime of the 
aligned lattice and the large- k regime of the anti-aligned 
lattice, the anomalous mode moves upward along the real 
axis as g m increases. In these cases, its collision with an- 
other mode yields a Hamiltonian-Hopf bifurcation and an 
oscillatory instability. In the large- k regime of the aligned 
lattice and the small- k regime of the anti-aligned lattice, 
the anomalous mode moves downward along the real axis 
as g m increases. When it hits the origin of the spectral 
plane of eigenfrequencies, one obtains an exponential in- 
stability associated with an imaginary eigenfrequency. 

We sweep over the parameter values g m G [0,1], 
g s G [0,0.95], and k G [0,15] and thereby confirm, for 
given g s and g m: that there is a critical wavenumber 
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FIG. 1: [Color online] Dark solitary- wave solutions of the 
stationary GP equation with an aligned nonlinear lattice 
(top panels) , their corresponding eigenfrequencies determined 
using the BdG equations (middle panels), and the spatio- 
temporal evolution of the solitary waves using simulations 
of the time-dependent GP equation (bottom panels) in an 
aligned nonlinear lattice. The left column shows an unstable 
dark solitary wave with a quartet of complex eigenfrequencies; 
the parameter values are g m = 1, g s = 0.5, k = 1, and \i — 10. 
The right column shows an unstable dark solitary wave with 
a pair of purely imaginary eigenfrequencies; the parameter 
values are g m = 1, g s = 0.5, k = 8, and \i — 10. In the top 
panels, we also show the harmonic trap (black dashed curves) 
and the inhomogeneity Ag(z) (red dash-dotted curves). We 
illustrate the spatio-temporal evolution of the BEC in the 
bottom panels; the color map indicates the value of \u(z,t)\ 2 . 




FIG. 2: [Color online] Dark solitary- wave solutions of the sta- 
tionary GP equation (top panels), their corresponding eigen- 
frequencies determined using the BdG equations (middle pan- 
els), and the spatio-temporal evolution of the solitary waves 
using simulations of the time-dependent GP equation (bottom 
panels) in an anti-aligned nonlinear lattice. The left column 
shows unstable solitary waves with a pair of purely imaginary 
eigenfrequencies; the parameter values are g m = 1, g s = 0.5, 
k = 1, and \i — 10. The right column shows an unstable 
dark solitary wave with a quartet of complex eigenfrequen- 
cies; the parameter values are g m = 1, g s = 0.5, k = 8, 
and fi = 10. In the top panels, we also show the harmonic 
trap (black dashed curves) and the inhomogeneity Ag(z) (red 
dash-dotted curves). We illustrate the spatio-temporal evolu- 
tion of the BEC in the bottom panels; the color map indicates 
the value of \u(z,t)\ 2 . 



k c for both aligned and anti-aligned nonlinear lattices. 
For aligned lattices, the anomalous mode increases as a 
function of g m away from g m = when k > fc c , and it 
decreases as a function of gm when k < k c . For anti- 
aligned lattices, however, the anomalous mode decreases 
as a function of g m away from g m = when k > fc c , and 
it increases as a function of g m when k < k c . We also 
trace the critical value k c between the small- k regime 
and the large- A: regime. We plot the dependence of k c 
on g s for fixed g m (top panels) and on g m for fixed g s 
(bottom panels) in Fig. for the aligned lattice and in 
Fig. [8] for the anti-aligned lattice. As one can see from 
the figures, the transition occurs in a narrow band near 
k = 6 for all values of g s and g m . Recall that the heal- 
ing length (which gives the length scale of the dark soli- 
tary wave) is £ = (87r7Voa) -1//2 , where No is the max- 
imum dimensional density of the 3D BEC. Specifically, 



No = ^o{Bhjj z /\q\)^ where q = 47rh 2 a/m = 2na 2 L q and 
no denotes the maximum density of the nondimensional 
solution. The wavelength of the lattice is 7ra z \/2/B /fc, so 
the ratio of the lattice wavelength to the healing length is 
ro = (7ra z /k)(167rNoa/ B) 1 / 2 . For example, with k — 6, 
the parameter value no ~ 4 yields ro ~ 2. As we discuss 
below, this corresponds to the critical region of param- 
eter space, when the width of the dark solitary wave is 
approximately two lattice wavelengths. Hence, when k 
is small, the variation of the scattering length occurs on 
a much larger scale than that of the solitary wave; when 
k is large, however, the variation occurs at a scale that 
is smaller than that of the solitary wave. Accordingly, a 
competition of length scales between the scale of the soli- 
tary wave and the scale of g(z) accounts for the existence 
of the two distinct regimes. 

For both aligned and anti-aligned lattices, the criti- 
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FIG. 3: [Color online] Variation of the eigenfrequency spec- 
trum as g m increases from for aligned lattices with g s = 0.5, 
k = 1, and \i — 10. The anomalous mode, which we plot using 
red circles, moves to larger frequencies as gm increases. From 
top to bottom, the values of gm are 0, 0.005, 0.008, and 0.015. 



FIG. 4: [Color online] Variation of the eigenfrequency spec- 
trum as g m increases from for aligned lattices with g s — 0.5, 
k = 8, and ji = 10. The anomalous mode, which we plot us- 
ing red circles, moves to lower frequencies (which eventually 
become imaginary) as g m increases. From top to bottom, the 
values of g m are 0, 0.005, 0.008, and 0.015. 



cal wavenumber k c exhibits a stronger dependence on g s 
than it does on Qm-) which is consistent with the fact that 
the depth of the wells depends increasingly sensitively 
on g s as g s — >• 1. In an aligned lattice with fixed g Sl 
the critical wavenumber k c decreases with g m , and this 
decrease becomes more dramatic as g s becomes smaller. 
By contrast, when g m is fixed, k c first decreases and then 
increases with g s ; this decrease becomes more dramatic 
as g m becomes smaller. This variation can be understood 
intuitively on the basis of length-scale competition. Sup- 
pose that ro ~ 2 for k = k c . An increase of g m then leads 
to a (maximal) density decrease and hence to a decrease 
of k c (because k c oc n^ 2 for ro ~ 2), and the opposite 
trend emerges from the increase of g s . In an anti-aligned 
lattice, these trends are reversed (as discussed above): 
with fixed g m , the critical wavenumber k c decreases with 
g Sl and this decrease becomes more dramatic as g m be- 
comes smaller. When g s is fixed, the BEC dynamics de- 
pends on whether g s is small (i.e., near 0) or large (i.e., 
near 1). When g s is small, k c stays almost constant when 



gm changes. When g s is large, however, k c increases with 

9m- 



IV. DYNAMICAL RESPONSE 

In Section [Till we examined the stability of dark soli- 
tary waves in BECs with a spatially periodic but time- 
independent nonlinearity coefficient g(z). We now con- 
sider BECs in which g is initially constant but then a 
Feshbach resonance is subsequently turned on to imple- 
ment a spatial dependence. This allows us to consider 
experimental situations in which the spatial dependence 
is turned on slowly (i.e., adiabatically) as well as ones in 
which it is turned on abruptly (i.e., nonadiabatically). 

We study the time-dependent GP equation (pQ) with 
g = g(z,t) and an initial wavefunction given by a sta- 
tionary solution to a GP equation with constant non- 
linearity coefficient g = go- The time-dependent and 
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FIG. 5: [Color online] Variation of the eigenfrequency spec- 
trum as gm increases from for anti-aligned lattices with 
g s — 0.5, k — 1, and \i — 10. The anomalous mode, which 
we plot using red circles, moves to lower frequencies (which 
eventually become imaginary) as gm increases. From top to 
bottom, the values of g m are 0, 0.007, 0.008, and 0.01. 



FIG. 6: [Color online] Variation of the eigenfrequency spec- 
trum as gm increases from for anti-aligned lattices with 
g s — 0.5, k — 8, and ji = 10. The anomalous mode, which 
we plot using red circles, moves to higher frequencies as g m 
increases. From top to bottom, the values of g m are 0, 0.002, 
0.005, and 0.007. 



space-dependent nonlinearity coefficient is given by 
g(z,t)=g + Ag(z)T(t), 



Tit) = \ 



1 + tanh 



t-t 



(3) 



where to denotes the time at which the function T(t) 
reaches the value 1/2. The function T(t) controlling the 
transition satisfies T(t) -)> as t -» and T(t) -)> +1 
as t — >• oo. For each transition time scale r, we choose 
to > 20t, which ensures that \T(t)\ is smaller than ma- 
chine precision (and hence approximately 0) at t = for 
our computations. This allows us to model the effect of 
the nonlinearity coefficient change from an initial value 
of go to a final value of go + Ag(z), where Ag(z) [see 
Eq. J2J)] represents the spatial dependence introduced by 
the Feshbach resonance. A small value of the parame- 
ter r corresponds to a situation in which the variation 
of the scattering length is abrupt, whereas a large value 
of r corresponds to a situation in which the Feshbach 



resonance is turned on gradually (i.e., adiabatically). 

In Fig. [9l we consider the transition time scales r = 0.1, 
r = 1, and r = 10 to illustrate fast, medium, and slow im- 
plementation of the spatial inhomogeneity in an aligned 
nonlinear lattice. Dark solitary waves appear to persist 
generically in our dynamical simulations. Nevertheless, 
unless the onset of the lattice is sufficiently adiabatic (i.e., 
for r = 10), additional excitations also arise. In partic- 
ular, when the nonlinear lattice is turned on sufficiently 
fast (see the panels corresponding to r = 0.1 and r = 1), 
we observe the emission of localized yet mobile gray soli- 
tary waves. Furthermore, for all three values of r, we 
observe oscillations of the dark solitary wave when we 
simulate long enough. For r = 1 and r = 10, the oscilla- 
tions start at t « 580 and t ~ 1.4 x 10 5 , respectively. This 
occurs at approximately t = to + 20r, which is when the 
nonlinear lattice has settled, within machine precision, 
to its final form. Dark solitary waves are stable if the 
nonlinearity coefficient is stationary and in its final state 
of g = go + Ag(z). 
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FIG. 7: [Color online] Dependence of the critical wavenumber 
k c on lattice parameter values for aligned nonlinear lattices. 
The top panel shows k c versus g s for multiple fixed values of 
gm, and the bottom panel shows k c versus g m for multiple 
fixed values of g s . The panels labeled "exponential instabil- 
ity" and "oscillatory instability" indicate regimes that exhibit 
those instabilities for the parameter values examined. 



The eventual mobility of the central solitary wave can 
be explained by its interaction with sound waves that are 
emitted during wave propagation [66]. The strength of 
the sound waves depends on the value of r. As r becomes 
larger, the system becomes more adiabatic, so the sound 
excitations become weaker and it takes longer for the 
above mechanism to destabilize the solitary wave. Given 
our collisionally inhomogeneous setup, eventual destabi- 
lization is inevitable. We plot the onset time of oscilla- 
tions T c = T c (t) in Fig. [TDJ By monitoring whether the 
oscillation has started before a given observation time 
T c , we can estimate a corresponding threshold time scale 
r c between nonadiabatic transitions and adiabatic tran- 
sitions. If the transition time is slower than r c , then 
no substantial oscillations occur before time t = T c . For 



FIG. 8: [Color online] Dependence of the critical wavenum- 
ber k c on lattice parameter values for anti-aligned nonlinear 
lattices. The top panel shows k c versus g s for multiple fixed 
values of gm, and the bottom panel shows k c versus gm for 
multiple fixed values of g s . 



T c = 500, T c = 1000, and T c = 5000, we obtain respective 
threshold time scales of r c « 0.6, r c « 2.8, and r c « 7.4. 
We test a series values of r and plot the associated T c 
values in Fig. [10l 



V. CONCLUSIONS 

In this paper, we have studied the structure, stabil- 
ity, and dynamics of dark solitary waves in collision- 
ally inhomogeneous Bose-Einstein condensates. We con- 
sidered spatially periodic scattering lengths (i.e., non- 
linear lattices) in a functional form suggested by re- 
cent experiments. Importantly, this family of nonlin- 
ear lattices can be tuned from a small-amplitude, ap- 
proximately sinusoidal structure to a periodic sequence 
of densely-spaced spikes. We demonstrated several in- 
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FIG. 9: [Color online] Dynamical response of dark solitary waves in BECs from introducing spatial heterogeneity into the 
nonlinearity coefficient by turning on a Feshbach resonance. From left to right, the columns correspond to simulations with 
transition time scales of r = 0.1 (nonadiabatic), r = 1 (medium), and r = 10 (adiabatic). The parameter values describing 
the shape of the nonlinear lattice are g m = 0.6, g s = 0.5, k = 0.5, and \i « 7.16. In the left panel of the first row, one can 
observe the emission of travelling gray solitary waves. We track one such wave: it reaches one of its left turning points (i.e., 
a point at which the wave's dip is at a value of z < such that \z\ is a local maximum) at t\ « 495, it then crosses z — at 
£2 ~ 518, and it reaches its next right turning point at £3 ~ 541. In the second row, we plot the corresponding spatial profiles 
of the solitary- wave solution of the GP equation at times £1 (left panel), £2 (center), and £3 (right). In the bottom panels, we 
circle the gray solitary that we track using a red dashed curve. The local minimum it represents is visible near z — —50 in the 
left column and near z — 50 in the right column. (It is contained within the dark solitary wave in the center column.) 
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FIG. 10: [Color online] Onset time T c for oscillations of a dark 
solitary wave versus the lattice transition time scale r c . 



teresting phenomena — including, for example, that dark 
solitary waves in aligned lattices and anti-aligned lat- 
tices exhibit different instability properties. These in- 
stability properties depend significantly on the modula- 
tion wavenumber of the nonlinear lattice, and the BEC 



dynamics exhibit considerable variation when the lattice 
wavelength is comparable to the size of the dark solitary 
waves. 

In the case of aligned lattices, we also examined the dy- 
namical response of solitary waves to adiabatic and nona- 
diabatic implementation of collisional inhomogeneities 
via a Feshbach resonance. When a Feshbach resonance is 
turned on nonadiabatically, nonlinear excitations can be 
emitted. Sound waves are emitted when a Feshbach res- 
onance is turned on (regardless of the transition speed), 
and these excitations eventually destablize the solitary 
wave. If the resonance is turned on sufficiently adiabati- 
cally, such that the sound waves are weak, then the BEC 
settles to an excited state (bearing a dark solitary wave) 
of the stationary nonlinear lattice initially before eventu- 
ally starting to oscillate after a very long time. 

It would be interesting to generalize our considerations 
to higher-dimensional lattices and to examine the effects 
that such lattices have on higher-dimensional excitations, 
such as vortices in quasi-2D BECs and vortex rings in 
3D condensates. Although some initial efforts have been 
made in that direction [64[, a systematic theory has yet 
to be developed, and thorough numerical investigations 
would also be beneficial. 
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